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Abstract 

We perform stability analysis of an accretion disc around a stellar mass black hole. 
The cooling of the accretion flow due to bremsstrahlung and synchrotron radia- 
tion processes are considered. The solutions with shock are perturbed at the shock 
location and the radial dependence of perturbation is obtained by solving the equa- 
tions numerically using fourth-order Runge-Kutta method. The perturbations are 
assumed to be infinitesimal, thereby allowing a linear analysis. The frequencies 
that solve the equations are found and the existence of unstable frequencies implies 
t^. , instability. 
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1. Introduction 



Most of the popular studies in accretion discs consider steady-state situations (Bondi 
1952; Shakura & Sunyaev 1973; Novikov & Thorne 1973; Narayan & Yi 1994; 
Paczyhski 1998) i.e. where things are not changing with time. Steady state ap- 
proximation is applicable only when the time scale for change of the physical phe- 
nomena is large. We will be better equipped to explain the transient and time- 
dependent phenomena by relaxing the assumption that the flow is steady in our 
analysis. The steady state properties of discs are mostly independent of viscosity 
and hence their observations do not lead to the quantification of the viscosity. But 
the time-dependent behaviour of the disc is more sensitive to the viscosity mecha- 
nism operating in the disc. Quantitative information of the viscosity can be derived 
by observations of the variable phenomena. An example of the time-dependent 
phenomena is the outburst of a dwarf novae. 

Analysing the stability of a fluid is one of the most difficult problems to handle in 
fluid dynamics. The monograph of Chandrasekhar (1961) uses variational principle 
which is an approximate method, to tackle the problem of hydrodynamic and hy- 
dromagnetic stability analytically. Langer, Chanmugam & Shaviv (1982, hereafter 
LCS82) consider column accretion in the case of a white dwarf in the presence of 
optically thin bremsstrahlung cooling. When the magnetic field is strong the matter 
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is accreted along the field lines onto the polar caps. The energy is released in form 
of extreme ultraviolet, soft X-ray and hard X-ray radiation. The time-dependent 
solution was obtained by a numerical simulation which uses finite difference scheme 
on an Eulerian grid. The post shock flow was found to be thermally unstable and 
the shock undergoes periodic oscillations. This oscillation was unexpected and in- 
teresting and the properties of the oscillation are dependent on the accretion rate 
and mass and radius of the white dwarf. The time scale of oscillation corresponds 
to post-shock cooling time scale. In Chanmugam, Langer & Shaviv (1985, hereafter 
CLS85) cooling due to cyclotron process also was included. The presence of cy- 
clotron cooling tends to stabilize the flow if the magnetic field is more than 2 x 10 7 
gauss. Balbus & Hawley (1991, hereafter BH91) consider magneto-hydrodynamics 
of accretion. The flow was found to be unstable to axisymmetric disturbances in 
the presence of a weak magnetic field. The instability is local and extremely pow- 
erful. The fluid motions create poloidal and toroidal magnetic field components. 
Chandrasekhar and Velikhov studied a related problem of Couette flow. This is a 
promising mechanism for generating viscosity but the derived viscosity coefficient 
value is not sufficient to explain the deduced accretion rates. 

There are different time scales in which the disc structure can vary. They are 
the viscous time scale t V i SC , dynamical time scale t^, time scale in which transverse 
hydrostatic equilibrium is established t z , thermal time scale t t h and cooling time 
scale t coo i. 

tvisc ~ — ~ — is the time scale for matter in the disc to diffuse due to viscous 
torques 

ts — — is the time scale of rotational motion 

t z = ^- is the time scale in which perturbations along z directions are smoothed 
out 

t th is the time scale for thermal perturbation to re-adjust 

t coo \ is the cooling time scale in which thermal energy is dissipated due to various 
cooling processes 

Oscillating shocks may have significance for accretion onto compact objects as it 
could potentially explain the observed X-ray variabilities, and in particular quasi- 
periodic oscillations (LCS82; Molteni, Sponholz & Chakrabarti 1996; Chakrabarti 
& Manickam 2000). We perform global stability analysis for an axisymmetric sys- 
tem and hence the perturbations take the form y 1 ( r ) e *( fcz + m< / > ^ w *) ) where /i(r) is an 
infinitesimal and it can be any one of the flow variables. 

2. Time-dependent infinitesimal perturbation from steady state 

We assume that the time-dependent solution can be expressed as a sum of time- 
independent solution and a time-dependent infinitesimal. Any arbitrary time-dependent 
infinitesimal can be expressed as sum of normal modes. We choose to express them 
as sines and cosines. The complex representation of the Fourier series is used to 
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study both the growth and decay of the modes. The variables are expressed as, 

p = p (r)+ Pl (r)e^ +m ^, 

v < i, = v < i >0 (r)+v <h (r)e i <- k * +m +- ut \ 
Vz=Vzi { r y(kz+ m ^t)^ 

and, 

p = Po {r) +p 1 (r)e i ( kz+m(t, -" t \ 

The unperturbed quantities have subscript and infinitesimals have subscript 1. The 
perturbations pi,v ri ,v ( f )1 ,v Zl ,pi are complex quantities. The wave number k takes 
real values, m is an integer (since + 2tt should be same as 0) and the frequency 
uj (= ujr + iujj) is a complex number. The presence of positive imaginary frequency 
(u>i) signals instability. 

3. Linear stability analysis 

The local stability analysis may miss an unstable mode which can be seen only in 
a global analysis. In the global analysis that is performed here, we assume that the 
perturbations are infinitesimals. So it is possible to neglect the higher order terms 
of perturbations other than the first order. This permits us to perform a linear 
analysis. For a non-linear analysis one has to rely on computer simulations. The 
surface r = r s where the Mach number reaches unity can be considered as a sound 
horizon, as the perturbations generated within that surface can be advected only 
inwards and it cannot propagate out. 

We use pseudo- Newtonian description of a Schwarzschild black hole with the help 
of a potential, prescribed by Paczyhski & Wiita (1980). This pseudo-potential pro- 
duces the positions of marginally stable and marginally bound orbits correctly and 
the efficiency approximately same as that of the effective potential of a Schwarzschild 
black hole. The equations describing the dynamics of the flow are the conservation 
equations of mass, momentum and energy. The variables velocity, density and pres- 
sure are denoted as v(iv, v^, v z ), p and p respectively. Cylindrical coordinate system 
(r, 0, z) is used and the time-dependent equations take the form, 

Continuity equation: 

dp 1 d 1 d d 

di + rd-r {rpVr) + rdfi" 3 ** + d~z ipVz) = °' M 
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The components of the Euler equation: 



_|_ v _|_ ^±^a + ^ ^ _ v _±_ _ _ 
<9t r <9r r 90 2 r p <9r <9r ' 

dvj, dv^ i^di^ dvf 1 1 dp 

~7tr + Vr A + id + Vz x + - VrV * = — *i> ( lc ) 

ar ar r <70 oz r pr o<p 



and, 



9i <9r r (90 z dz pdz' 



Energy equation: 



de de de de , . 

pVr dr + pV<t> ~rd4> + pVz ~dz + brems + sync = ~ P di' 

where, 

g = _ CM. is 

the Paczyhski-Wiita pseudo-Newtonian potential, G is the 

gravitational constant, M is the mass of the black hole, r g = ^ is the Schwarzschild 
radius and c is the velocity of light. The specific energy(e) of the flow is \{v^ + + 
v z) + pt^T 9 m ^ ne ex P ress i° n f° r bremsstrahlung (Lang 1980) and synchrotron 
(Shapiro & Teukolsky 1983) cooling in an electron-proton plasma are, 

A brems = 1.43 x 1 - 27 4t 1/2 S/, 

mi 



and, 

3 c m e c 7(7 — 1) m e cr 7 m p 

respectively, where m p is the mass of proton, T is the temperature, gf is the Gaunt 
factor, e is the electron charge, m e is the electron mass, 7 is the adiabatic index 
and p = 0.5 for purely hydrogen gas. Using the ideal gas equation, we get for the 
temperature T, 

= p pm p 
p k 

where, k is the Boltzmann constant. The polytropic relation p = Kp 1 is used to 
obtain the sound speed a = = 

Frequencies that exist in the system are determined by a 'dispersion relation'. 
Here in this case we do not have an explicit equation for the dispersion relation 
(See BH91 for an example of an explicit dispersion relation). The set of equations 
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which are obtained for infinitesimals can be considered as 'dispersion relation'. The 
frequencies that solve this set of equations are the eigen frequencies which exist in 
the physical system. The variables of the perturbed flow are expressed as a sum of 
unperturbed value (which is the steady state solution) and an infinitesimal. These 
variables are introduced in the flow equations and simplified ignoring higher order 
terms. The equations are split into real and imaginary parts to obtain following ten 
equations involving ten unknowns. 

T A V' TR = ~UJ R T 2 - W/Tg - V ro v' ro V rR - VroV^Vfa - UroWto + ^iK.^K) ~ 

(2) 

T 4 v' rj = -uj r T 3 -uj{T 2 -v ro v' ro v ri -v ro v'^ v (j}! -v ro v'^ + ^ ( v ro ^p' -2v ro p lPo 3j ) + 
^v ro f/ - v ri Ti + >T 3 - 4^(f )^^PoTo\{f + %)- ^% - ^r (-Ts + 

2Li 



(3) 

T 5 + PrVvo + Pov' rR = 0, (4) 

T 7 + p'jV ro + p v' ri = 0, (5) 



T 6 + VroV' = - — {p' R - —Po), (6) 

Po Po 

Ts + v ro v' ri ^--(p' I -^p' ), (7) 
Po Po 



i i v 4>o 1 1 1 1 
+ VfaUR + VroVj + VrnVfo V^m + -V rQ V <j>R + -V rR V^ = P/771, (8) 

r r r p r 



- f^c^fl + v ro v^ + Urj^,, + —v<t> R m + -v ro v <t>I + -v^v^ = - — -p R m, (9) 

v Zr uji + v 2/ cjr + Uro^L -v Zl m = —pik, (10) 

r Po 



and, 



u^wj - v Zr uj r + u^u^ + —v ZR m = — -p R k, (11) 

r Po 
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where, 



rp , , / 7 9 Pq , 

Ti = v ro v ro + + 7— y ^(-) + 5 , 



T 2 = u u + + — - ( p + Pl ), 

17 - IJPo Po 

T 3 = V ro V rR + VfoVfa + ~ r— ( Po + PR), 

(7 - IJPo Po 



2 7 Po 7 2 

7 V 

ro 7-lpo 7-1 ro ' 



1 / 1 1 , 1 

X 5 = PR^/ + P/^i? + -Pi?^r + Pfl^ro ~~ -PoV^m - -piTTlVfo - p V Zl k + -p Q V rR + p Q V rR , 

T 6 = w ri{ cj/ + t^OJfl + Wr fl ^ ro y v r 1 , m — v 4>r-. 



1 / 1 1 , 1 

T 7 = p/^/ - p R uj R + -p/^r + p/^ ro + -pov^m + -pR-rnvfo + p v ZR k + -pov ri + p v ri , 



and, 



rr , / , v 4>o 2v <f>o 
T 8 = v ri uj! - v rR u R + v ri v ro H v rR m v^. 



r r 



Here, ' denotes the derivative with respect to r. Fourth-order Runge-Kutta method 
(Press et al. 1992) is used to solve the equations (2-11). The perturbations are 
introduced at the shock location and the equations are integrated to obtain the 
functions, pr, pi, v rR , v ri , v^ R , v^ n v ZR , v Zl , pr and pj. The frequencies are changed 
over a grid. 



4. Eigen frequencies 

The 'dispersion relation' gives the frequencies that exist in the system. All the 
quantities oscillate with that frequency but with different amplitudes. As the shock 
location oscillates, the size of the post-shock region oscillates thereby intercepting 
different amounts of soft X-ray flux, emitted by the pre-shock disc of a stellar mass 
black hole. The fractional change in size of the post-shock region is large compared 
to that of the pre-shock disc. Hence QPOs are expected mainly in the hard X-ray 
flux. For the case of super-massive black hole, the post-shock region radiates at a 
lower frequency. 
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We choose the mass of the stellar black hole to be 14M Q and the value of 7 
equal to 4/3, as distinct X-ray variability of GRS 1915+105 is interpreted as due to 
instabilities in radiation pressure dominated disc (Greiner, Cuby and McCaughrean 
2001). The instability is decided by the presence of positive imaginary part of 
the perturbation frequency u>j. In Manickam (2004, hereafter Paper I), the flow 
topologies and the parameter dependence on the possibility of the formation of 
Rankine-Hugoniot shock was studied. It was mentioned that among the different 
branches that are possible, for a given accretion rate and specific angular momentum 
(A), as a valid shock solution, stability and boundary conditions should decide the 
uniqueness of the solution. Here we consider the case, where accretion rate is equal 
to unity in Eddington units and specific angular momentum is 1.8 in the units 
of 2GM/c. The stability is analysed for bremsstrahlung cooling case and when 
both bremsstrahlung and synchrotron cooling are present. The efficiency factor for 
synchrotron cooling process is chosen as 10~ 7 and 10~ 6 , while that of bremsstrahlung 
cooling process is chosen as unity. 

The branches of flow topologies that are considered for stability analysis are 
shown in Figs, l(a-c). The stability of different branches when perturbations are 
introduced at outer shock location are analysed. The branches which form a shock 
are shown in r shk vs r cout plot of Fig. 2. We choose the supersonic branches cor- 
responding to the outer critical point location at 23r g , 25r g and 27r g for stability 
analysis. The frequencies that exist in the system are obtained by numerically in- 
tegrating the ten equations, (2-11), using fourth-order Runge-Kutta method (Press 
et al. 1992). The following is the 'flow chart' which is implemented in the code. 

(i) the mass of the black hole is chosen to be 14M Q , and 7 as 4/3 

(ii) set the efficiency factor of bremsstrahlung cooling (zefac) and synchrotron 
cooling (synfac) 

(iii) obtain the flow topology (as described in Paper I) and find out the branches 
that will form a shock 

(iv) choose an outer supersonic branch and an inner subsonic branch that satisfies 
the shock conditions 

(v) choose the real and imaginary frequencies of perturbation, ujr and ui 

(vi) introduce perturbation at the shock location 

(v) integrate using fourth-order Runge-Kutta method to obtain the amplitude of 
the perturbation in the subsonic branch (supersonic branches are not considered 
for stability analysis as the perturbations would get advected with the flow in 
infall time scale) 

(vi) frequencies are changed over a grid 

(vii) if a positive imaginary frequency exists, then the accretion disc is unstable 

(viii) if perturbation amplitude blows up for even one of the flow variables, then 
the flow is not amenable to a linear analysis 
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5. Results 

The Fig. 3(a-b), Fig. 4(a-b) and Fig. 5(a-b) show the frequencies that exist in the 
system. These are plots of (uj/2n) vs (ujr/I'k) . Hence values in the axis are in 
hertz. It is of the order of 10 Hz which is same as that observed in galactic black 
hole candidate GRS 1915+105. For the case of synchrotron efficiency factor of 10~ 7 , 
there is not much effect in the frequency modes that exist in the system, due to 
addition of synchrotron cooling process. 

But when the efficiency factor is increased to 10~ 6 , the flow topology changes 
drastically (Fig. lc) and the flow is free of shock. This is in line with the observa- 
tion of CLS85, where shock disappears if the magnetic field is above a certain value. 
Stability analysis with more realistic models could explain the observed QPOs bet- 
ter. For instance, the effect of magnetic field on the dynamics of the flow is not 
considered. Magnetic field seems to be present in GRS 1915+105 and it plays 
an important role (Nandi, Chakrabarti, Vadawale & Rao 2001; Vadawale, Rao & 
Chakrabarti 2001). 

6. Discussion and Conclusions 

Not all the frequency ranges fall within the linear analysis domain, and the failure 
of linear analysis does not imply instability or stability. Like the accreting fluid has 
infinite degrees of freedom, there are many cases which can be considered for linear 
analysis. Here only a few representative cases are illustrated. These results could 
be useful to test a time-dependent numerical code which could handle non-linear 
perturbation analysis. 
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Fig. la: The branches of the flow topology considered for stability analysis. Su- 
personic branches have critical point locations at 23r g , 25r g and 27r g . Efficiency of 
bremsstrahlung cooling is unity and there is no synchrotron cooling. Mass of stellar 
black hole is 14M and 7 = 4/3. Accretion rate is 1.0 in Eddington units and A=1.8 
in units of 2GM/c. 
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Fig. lb: Same as Fig. la but synchrotron cooling of efficiency factor 10 7 is included. 
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Fig. lc: Same as Fig. la but synchrotron cooling of efficiency factor 10 6 is included. 
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Fig. 2: Outer critical point location (r cout ) of the supersonic branches that forms 
shock at location r s hk- The left set of points correspond to bremsstrahlung plus 1CT 7 
efficient synchrotron cooling process and right set is for bremsstrahlung cooling. The 
parameters are same as that of Figs, l(a-b). 
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Fig. 3(a-b): This shows the frequency modes that exist in the accretion disc for 
outer supersonic branch of critical point location at 23r g . The positive imaginary 
frequencies are the unstable modes. Bremsstrahlung cooling case is shown in (a) 
and bremsstrahlung plus 10~ 7 efficient synchrotron cooling is shown in (b). 
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Fig. 4(a-b): Same as Fig. 3(a-b) but outer supersonic branch has critical point 
location at 25r g . 
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Fig. 5(a-b): Same as Fig. 3(a-b) but outer supersonic branch has critical point 
location at 27r g . 
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